/*
 * Author  : Peter Mosbach & Keith Vetter
 * Date    : August 2003
 * Desc    : Math functions for equation parser.
 *           The inverse error algorithm was originally
 *           a piece of FORTRAN 90 code that we found on the
 *           internet.
 */

#include <math.h>

#include "eqparse_protos.h"


double add(double x,double y)
{
	return(x+y);
}

double dif(double x,double y)
{
	return(x-y);
}

double mul(double x,double y)
{
	return(x*y);
}

double divide(double x,double y)
{
	if (y!=0) {
		return(x/y);
       } else {
              return(0) ;
       }
}

double mod(double x,double y)
{
	return(fmod(x,y));
}

/* Inverse error function */
double inverf(double p)
{

        double a[6]  = {  .1400216916161353e+03, -.7204275515686407e+03, 
                          .1296708621660511e+04, -.9697932901514031e+03,
                          .2762427049269425e+03, -.2012940180552054e+02 } ; 

        double b[6] = { .1291046303114685e+03, -.7312308064260973e+03,  
                        .1494970492915789e+04, -.1337793793683419e+04, 
                        .5033747142783567e+03, -.6220205554529216e+02 } ;
        
        double a1[7]  = { -.1690478046781745e+00, .3524374318100228e+0,  
                          -.2698143370550352e+02, .9340783041018743e+02,  
                          -.1455364428646732e+03, .8805852004723659e+02,  
                          -.1349018591231947e+02 };

        double b1[7] =  { -.1203221171313429e+00, .2684812231556632e+00,
                          -.2242485268704865e+02, .8723495028643494e+02, 
                          -.1604352408444319e+03, .1259117982101525e+03, 
                          -.3184861786248824e+02 } ;

        double a2[9] = { .3100808562552958e-04, .4097487603011940e-02,
                         .1214902662897276e+00, .1109167694639028e+01, 
                         .3228379855663924e+01, .2881691815651599e+01, 
                         .2047972087262996e+01, .8545922081972148e+00,
                         .3551095884622383e-02 };

        double b2[8] = { .3100809298564522e-04, .4097528678663915e-02,
                         .1215907800748757e+00, .1118627167631696e+01,
                         .3432363984305290e+01, .4140284677116202e+01,  
                         .4119797271272204e+01, .2162961962641435e+01 };

        double a3[9] = { .3205405422062050e-08, .1899479322632128e-05,
                         .2814223189858532e-03, .1370504879067817e-01,
                         .2268143542005976e+00, .1098421959892340e+01,
                         .6791143397056208e+00, -.834334189167721e+00,
                         .3421951267240343e+00 } ;

        double b3[6] = { .3205405053282398e-08, .1899480592260143e-05,
                         .2814349691098940e-03, .1371092249602266e-01,
                         .2275172815174473e+00, .1125348514036959e+01 };

        double  c = .5625  ;
        double c1 = .87890625 ;
        double c2 = -.2302585092994046e+03 ;
        double  eps, q, s, t, v, v1 ;

        if ( p<0.5 ) {
                q=p;
        } else {
                q=1.0-p;
        }

        if (p >= 0.0 && q > 0.0) {

                /*eps = MAX(EPSILON(1.0), 1.e-15) F90 code*/
                eps = 1.0e-15 ;
                t = 0.5 + (0.5-(p+q)) ;

                if (fabs(t) > 3.0*eps) {
                        return( -2.0 ) ;
                }
        
        
                if (p <= 0.75) {
                        v = p * p - c ;
                        t = p * (((((a[5]*v + a[4])*v + a[3])*v + 
                                     a[2])*v + a[1])*v + a[0]) ;
                        s = (((((v + b[5])*v + b[4])*v + b[3])*v + 
                                 b[2])*v + b[1]) * v + b[0] ;
                } else {
                        
                        if (p <= 0.9375) {
                                v = p * p - c1 ;
                                t = p * ((((((a1[6]*v + a1[5])*v + a1[4])*v + 
                                    a1[3])*v + a1[2])*v + a1[1])*v + a1[0]) ;
                                s = ((((((v + b1[6])*v + b1[5])*v + b1[4])*v + 
                                    b1[3])*v + b1[2])*v + b1[1]) * v + b1[0] ;
                        } else {
        
                                v1 = log(q) ;
                                v = 1.0 / sqrt(-v1) ;
                                if (v1 >= c2) {
                                        t = (((((((a2[8]*v + a2[7])*v + 
                                             a2[6])*v + 
                                             a2[5])*v + a2[4])*v + a2[3])*v + 
                                             a2[2])*v + a2[1]) * v + a2[0] ;
        
                                        s = v * ((((((((v + b2[7])*v + 
                                            b2[6])*v + 
                                            b2[5])*v + b2[4])*v + b2[3])*v + 
                                            b2[2])*v + b2[1])*v + b2[0]) ;
                                } else {
                                        t = (((((((a3[8]*v + a3[7])*v + 
                                            a3[6])*v + 
                                            a3[5])*v + a3[4])*v + a3[3])*v + 
                                            a3[2])*v + a3[1]) * v + a3[0] ;
        
                                        s = v * ((((((v + b3[5])*v + 
                                                 b3[4])*v +
                                                 b3[3])*v + b3[2])*v + 
                                                 b3[1])*v + 
                                                 b3[0]) ;
                                }
                        }
                }
        
                return( t/s ) ;
        }

        return( -1 ) ;
}
